Raster data in R
This document is a ‘cheatsheet’ for working with raster datasets in R and introduces the terra package which is the preferred option for handling such data currently. We first load some libraries (including some vector libraries too).
library(terra)
library(sf)
library(tmap)
library(dplyr)
library(tidyr)
library(ggplot2)
tmap_mode("plot")An example dataset to start
Reading a dataset is pretty simple, especially if it is a simple dataset. I have put a bunch of GeoTIFFs in this zip file, which you should download and unzip to a folder on your computer. Then read the digital elevation model dataset dem.tif into a variable with the rast() function, like this:
dem <- rast("dem.tif")Confusingly this doesn’t tell you anything about what it just did. So you have to ask:
dem## class : SpatRaster
## dimensions : 242, 240, 1 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 1735056, 1759056, 5419610, 5443810 (xmin, xmax, ymin, ymax)
## coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization)
## source : dem.tif
## name : dem
## min value : 0
## max value : 518.799This tells us a few things about the dataset, including the dimension (in numbers of cells), the resolution (in units of the projection), the extent, and projection . It also tells us the minimum and maximum values in the data.
As usual in R we can get a summary of the statistics of the dataset:
summary(dem)## dem
## Min. : 0.00
## 1st Qu.: 82.07
## Median :149.17
## Mean :156.96
## 3rd Qu.:219.22
## Max. :518.80
## NA's :23790What are all those NA values? Well… best bet is to plot the data to make sense of it
plot(dem)
This shows us right away that a lot of the rectangular area of the raster is taken up by sea, and since it is an elevation model, we might expect the sea areas to be recorded as NA values (which they are).
tmap can consume rasters too!
tmap knows about raster data and can map them, quite happily. The same options such as the number of classes (n), the style of classification (e.g. quantile), and the colour palette to us in the mapping also work for rasters. Here, because you might want to use it, I’ve also shown how to use R’s built-in terrain colour palette.
tm_shape(dem) +
tm_raster(n = 10, palette = terrain.colors(10)) +
tm_legend(legend.outside = TRUE)## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
tmap can also map raster layers in combination with non-raster layers. For example, read in the nz.gpkg file also found in the provided zip file, and add it as an outline around the raster:
welly <- st_read("welly.gpkg")## Reading layer `welly' from data source
## `/home/osullid3/Documents/teaching/Geog315/labs/scratch/rasters/welly.gpkg'
## using driver `GPKG'
## Simple feature collection with 79 features and 64 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1735078 ymin: 5419585 xmax: 1759042 ymax: 5443764
## Projected CRS: NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization)tm_shape(dem) +
tm_raster() +
tm_shape(welly) +
tm_borders(col = "blue")## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
But wait, there’s more than one layer!
There are a bunch of raster data sets in this folder. We can read them into a raster stack by making a directory listing of them and supplying that to the rast function
layer_sources <- dir(pattern = "*.tif")
layers <- rast(layer_sources)Again, you have to ask to get more information
layers## class : SpatRaster
## dimensions : 242, 240, 11 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 1735056, 1759056, 5419610, 5443810 (xmin, xmax, ymin, ymax)
## coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization)
## sources : age.tif
## deficit.tif
## dem.tif
## ... and 8 more source(s)
## names : age, deficit, dem, mas, mat, r2pet, ...
## min values : 0.0000, 0.0000, 0.0000, 138.0000, 100.0000, 24.0000, ...
## max values : 2.00000, 133.02213, 518.79901, 143.00000, 131.00000, 41.00000, ...If we want to know the names of all the layers, we can use the names() function
names(layers)## [1] "age" "deficit" "dem" "mas" "mat" "r2pet" "rain"
## [8] "slope" "sseas" "tseas" "vpd"And plot will make several plots of all the layers (this might be a bit slow).
plot(layers)
tmap can also deal with multiple layers, although it will apply the same classification to all of them unless you override this behaviour:
tm_shape(layers) +
tm_raster() + #title = names(layers)) # this fix no longer needed
tm_layout(legend.position = c("RIGHT", "BOTTOM")) +
tm_facets(free.scales = TRUE) # different scales on each layer
There might be some complaining about fitting the legends in on this plot, or about projections (this is to do with recent changes in how projections are handled by the proj tools see this post but you can ignore it for now).
To examine a single layer in the stack, just use the $ notation:
tm_shape(layers$mat) +
tm_raster(palette = "-RdBu", n = 9, title = "Mean annual temp (x10)")## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
Layers in the raster stack
The layers included now in the layers dataset are:
The layers in the stack were assembled from LENZ data developed by Manaaki Whenua Landcare Research, and available here. Their interpretations are roughly as follows (this information from the disdat package).
| Layer name | Explanation | Units |
|---|---|---|
| age | Soil age | 3 classes (0 to 2): <2000, 2000-postglacial (app. 30,000), and pre-glacial |
| deficit | Annual water deficit | mm (higher is greater deficit) |
| dem | Elevation | metres |
| mas | Mean annual solar radiation | MJ/m2/day |
| mat | Mean annual temperature | degrees C * 10 |
| r2pet | Average monthly ratio of rainfall and potential evapotranspiration | (ratio) |
| rain | annual precipitation | mm |
| slope | Slope | degrees |
| sseas | Solar radiation seasonality | dimensionless |
| tseas | Temperature seasonality | degrees C |
| vpd | Mean October vapor pressure deficit at 9 AM | kPa |
Raster calculation
terra will also allow you to do operations between layers, to perform raster calculations. For example, we can easily do this (nonsensical) calculation
layers$mas + layers$mat - layers$dem * layers$deficit## class : SpatRaster
## dimensions : 242, 240, 1 (nrow, ncol, nlyr)
## resolution : 100, 100 (x, y)
## extent : 1735056, 1759056, 5419610, 5443810 (xmin, xmax, ymin, ymax)
## coord. ref. : NZGD2000 / New Zealand Transverse Mercator 2000 (with axis order normalized for visualization)
## source : memory
## name : mas
## min value : -14228.81
## max value : 274This makes a new layer by mathematically combining the values in the original layers. When combined in this way raster layers must be aligned with one another and in the same projected coordinate system, and the values at the same locations in each layer contribute to the result at the corresponding location in the result.
This means that we can use raster calculations of this type to perform overlay as described in the next section.
Raster overlay
If we recode each layer by whether or not it is greater than some threshold, then we can logically combine them to find regions that satisfy a desired combination of criteria. Say we wanted places below 250m elevation with slope under 5 degrees:
below_250m <- layers$dem < 250
flat <- layers$slope < 5
result <- below_250m & flatAs usual make some maps
m1 <- tm_shape(below_250m) +
tm_raster(palette = "Set1", style = "cat", title = "dem<250m") +
tm_shape(welly) +
tm_borders(col = "white", lwd = 0.2)
m2 <- tm_shape(flat) +
tm_raster(palette = "Set1", style = "cat", title = "slope<5deg") +
tm_shape(welly) +
tm_borders(col = "white", lwd = 0.2)
m3 <- tm_shape(result) +
tm_raster(palette = "Set1", style = "cat") +
tm_shape(welly) +
tm_borders(col = "white", lwd = 0.2)
tmap_arrange(m1, m2, m3, nrow = 1)
There is nothing to stop us performing all these operations in a single step, once we are used to how it works:
r <- (layers$dem < 250) & (layers$slope < 5)This makes it relatively easy to flexibly perform complex suitability evaluations using terra raster calculation operations.
Focal calculations
The simple calculations in the previous sections are known as local because they apply locally across the values at the same locations in each layer. Two other kinds of calculation are possible focal and zonal.
We can apply functions in focal windows within a single layer and perform calculations based on the values within that window centred on each location across a raster. For example, we might want to determine the max value of raster cells in a 3x3 window.
max_elev <- layers$dem %>%
focal(fun = max)
plot(max_elev)
A 3x3 window is the default. If you want a larger window specify the size with the w parameter (which must be an odd number)
max_elev <- layers$dem %>%
focal(w = 11, fun = max)
plot(max_elev)
Notice how a larger window produces a more ‘smoothed’ or generalised result.
Again, we can flexibly combine such operations. For example a measure of relief (or roughness) is the local difference between the maximum and minimum values of that surface
relief <- focal(layers$dem, fun = max) - focal(layers$dem, fun = min)
names(relief) <- "relief"
tm_shape(relief) +
tm_raster(style = "cont", palette = "Reds", alpha = 0.7)## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
Zonal calculations
You can also perform calculations based on a set of zones. If you want the results in a raster, then both inputs must be rasters. That means we have to do quite a lot of fiddling around with the data. First we make a zones dataset from some polygons, converting first to terra‘s SpatVector format, then using rasterize to select one attribute from the zones as an identifier for zones. The raster we provide to the rasterize function is a ’template’ for the number of cells and their resolution.
welly_zones <- welly %>%
as("SpatVector") %>%
rasterize(relief, field = "ID")
tm_shape(welly_zones) +
tm_raster(n = 30, palette = "Set1", legend.show = FALSE) +
tm_shape(welly) +
tm_borders()## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
We can’t easily map the zones, because there are a large number of unique values.
Once we have zones we can perform a zonal calculation, which will calculate the result of applying a function across all the raster cells in each zone. We may need to change the NA values in the input raster layer to zeros (or some other defined value) so that any missing NA values in a zone don’t prevent calculation of a result for the zone they are in.
zonal_relief <- relief %>%
replace_na(0) %>%
zonal(welly_zones, fun = mean, as.raster = TRUE) %>%
mask(relief) ## The result of this operation is a raster where the value in each cell is the mean relief of all the cells in the zone it is in.
tm_shape(zonal_relief) +
tm_raster(n = 10, palette = "RdYlBu", title = "Zonal mean relief") +
tm_shape(welly) +
tm_borders()## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
In many situations, it is more appropriate to output the results from zonal calculations as a table of values which can be attached to the original vector data.
zonal_relief_table <- relief %>%
replace_na(0) %>%
zonal(welly_zones, fun = mean)
welly_relief <- welly %>%
merge(zonal_relief_table)
tm_shape(welly_relief) +
tm_polygons(col = "relief", style = "cont")
Surface analysis
terra provides an array of commonly used surface analysis methods. These are a better option than figuring out how to calculate them yourself using focal functions.
For example, we can calculate slopes from an elevation model using the terrain function
slope <- terrain(dem)
tm_shape(slope) +
tm_raster()## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
The terrain function offers several options, which you can find out more about in the help using ?terrain. A common use is for calculating hillshade, which can be done like the example below
aspect <- dem %>%
terrain("aspect", unit = "radians")
slope <- dem %>%
terrain("slope", unit = "radians")
hillshade <- shade(slope, aspect, angle = 315)
tm_shape(hillshade) +
tm_raster(palette = "Greys", style = "cont", alpha = 0.65, legend.show = FALSE)## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition
Combining raster and vector operations
In the examples for zonal calculations, it was necessary to combine raster and vector data and translate between them.
This is a common requirement, and can sometimes be challenging. Here a few important examples are shown. The key to success is realising that terra has its own vector data format SpatVector and terra functions want vector data to be in this format in order to work. To convert sf vector data to SpatVector use as("SpatVector") and to convert from SpatVector to sf use st_as_sf().
Extract values to points
For example, say we have a set of points, which we want values from a raster layer attached to.
First load the points. Here we have some sampled plants. They have to be projected to match the target raster, and we also restrict them only to the study area using st_filter().
plants <- st_read("tradescantia.gpkg") %>%
st_transform(crs(layers)) %>%
st_filter(welly)## Reading layer `tradescantia' from data source
## `/home/osullid3/Documents/teaching/Geog315/labs/scratch/rasters/tradescantia.gpkg'
## using driver `GPKG'
## Simple feature collection with 1635 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -176.5629 ymin: -46.9931 xmax: 178.4418 ymax: -34.7562
## Geodetic CRS: WGS 84Extract the values of a particular raster with this point set, we have to convert it to SpatVector points, then run terra::extract (we have to specify terra:: because there is another function called extract in tidyr).
mean_t <- layers$mat %>%
terra::extract(as(plants, "SpatVector"))
as_tibble(mean_t)## # A tibble: 120 × 2
## ID mat
## <dbl> <dbl>
## 1 1 125.
## 2 2 119
## 3 3 126.
## 4 4 114.
## 5 5 117
## 6 6 130
## 7 7 121
## 8 8 123.
## 9 9 125
## 10 10 120.
## # … with 110 more rowsThe result is just a table. If we want to map this, we have to add it back to the original spatial data:
plants <- plants %>%
mutate(mean_t = mean_t$mat)
tm_shape(plants) +
tm_dots(col = "mean_t")
We can apply terra::extract to a raster stack to get a multi-attribute table all in one go:
layers %>% terra::extract(as(plants, "SpatVector")) ## ID age deficit dem mas mat r2pet rain slope sseas
## 1 1 2 62.28255 81.850563 142 125.0223 29.00000 1181.559 24.977749 44
## 2 2 2 26.64497 170.385193 141 119.0000 35.00000 1206.175 7.644967 44
## 3 3 2 48.87741 74.046143 141 125.9484 32.05157 1193.613 23.967861 44
## 4 4 2 14.00000 270.647339 139 113.9835 39.00000 1171.995 17.677706 44
## 5 5 2 25.39118 200.619965 141 117.0000 36.00000 1153.038 6.008887 44
## 6 6 2 98.96744 8.207642 143 130.0000 26.00000 1115.785 1.000000 45
## 7 7 2 49.60237 141.479630 142 121.0000 32.00000 1134.068 14.609330 44
## 8 8 2 41.35058 115.502556 140 122.9810 35.00000 1158.736 15.000000 44
## 9 9 2 45.62231 83.624641 141 125.0000 33.00000 1190.220 13.034999 44
## 10 10 2 26.95165 167.459930 140 119.9517 36.00000 1197.365 15.633889 44
## 11 11 2 44.27446 92.351151 141 124.9688 33.00000 1193.510 24.991325 44
## 12 12 2 45.64238 75.081467 141 125.0000 34.00000 1183.438 9.397834 44
## 13 13 2 35.00000 131.923416 141 122.0000 34.00000 1201.759 2.625707 44
## 14 14 2 109.00000 2.933569 143 130.0000 25.00000 1137.889 5.926316 45
## 15 15 2 68.56382 50.979637 142 127.0000 29.00000 1148.681 12.436181 44
## 16 16 2 66.77875 78.654251 142 125.0000 29.00000 1131.083 15.435949 44
## 17 17 2 14.00681 244.173737 140 115.0068 39.00000 1160.724 17.981869 44
## 18 18 2 30.37409 154.893311 141 120.0000 34.00000 1224.243 14.041884 44
## 19 19 2 52.31923 121.309883 142 122.6391 31.00000 1159.858 23.655590 44
## 20 20 2 13.62304 230.870667 140 115.0265 38.00000 1183.262 13.911610 44
## 21 21 2 42.73308 105.846130 141 123.0000 34.00000 1163.599 12.867723 44
## 22 22 2 31.01176 138.779724 141 121.0000 36.00000 1171.339 5.000000 44
## 23 23 2 53.10150 45.262615 142 127.0508 32.00000 1192.101 21.796999 44
## 24 24 2 32.07914 150.639267 141 121.0000 34.00000 1234.301 5.029380 44
## 25 25 2 37.02601 114.676880 141 123.0000 35.00000 1184.033 11.641457 44
## 26 26 2 41.03801 90.667397 141 124.0000 34.00000 1181.224 9.061888 44
## 27 27 2 42.03058 123.655266 141 122.0306 34.00000 1163.380 7.030577 44
## 28 28 2 44.29372 83.716370 141 125.0000 34.00000 1182.644 15.311479 44
## 29 29 2 77.61778 8.234794 142 129.0000 29.00000 1154.349 13.351653 44
## 30 30 2 42.03058 123.655266 141 122.0306 34.00000 1163.380 7.030577 44
## 31 31 2 65.62665 83.118309 141 123.9973 32.00000 1156.229 15.753992 44
## 32 32 2 79.00000 8.101830 142 130.0000 29.00000 1152.648 1.000000 44
## 33 33 2 53.00000 51.516697 141 127.0000 32.00000 1190.664 9.076740 44
## 34 34 2 49.05710 111.361656 141 123.0000 32.00000 1192.307 8.326986 44
## 35 35 2 70.97032 65.887962 142 126.0000 29.00000 1128.287 7.624279 44
## 36 36 2 73.60202 28.743425 142 128.9852 29.00000 1150.323 8.360874 44
## 37 37 2 34.58276 132.193329 141 121.9512 35.00000 1196.248 14.582765 44
## 38 38 2 27.04408 161.311722 141 119.3645 35.03992 1178.057 27.404453 44
## 39 39 NaN NaN NaN NaN NaN NaN 1173.061 NaN NaN
## 40 40 2 60.37127 64.083328 142 126.0000 31.00000 1154.409 8.419346 44
## 41 41 2 53.52137 55.878609 142 126.9808 32.00000 1193.438 20.835752 44
## 42 42 2 83.57523 52.906380 142 127.0000 27.00000 1119.425 12.982502 45
## 43 43 2 74.15147 88.796715 142 124.9903 28.00000 1108.978 14.013801 45
## 44 44 2 21.00751 184.083023 140 118.0000 37.00000 1161.571 5.980017 44
## 45 45 2 14.00000 221.800766 140 116.0000 39.00000 1169.772 10.629755 44
## 46 46 2 41.76005 106.200203 141 123.9693 33.00000 1189.269 20.429573 44
## 47 47 2 57.42181 61.602367 142 127.0000 32.00000 1168.331 10.370629 44
## 48 48 2 52.00000 121.581322 142 122.0000 32.00000 1137.967 6.771970 44
## 49 49 2 83.73991 30.918926 140 128.0462 30.00000 1192.443 18.597193 44
## 50 50 2 28.01720 179.246902 141 118.0000 36.00000 1153.922 10.371836 44
## 51 51 2 26.68111 159.442719 141 120.0000 36.00000 1182.651 14.334340 44
## 52 52 2 47.61535 97.040489 141 124.0000 33.00000 1162.150 8.384654 44
## 53 53 2 29.00000 156.824371 141 120.0000 35.00000 1208.444 4.940795 44
## 54 54 2 33.77411 166.106613 141 119.0000 35.00000 1149.401 19.612944 44
## 55 55 2 48.94762 84.229973 141 124.9676 33.00000 1166.231 10.333567 44
## 56 56 2 45.33219 90.564613 141 124.9517 33.00000 1187.568 18.712646 44
## 57 57 2 36.40050 181.974335 141 118.0074 34.00000 1141.459 14.385720 44
## 58 58 2 45.24791 18.098948 142 128.3530 32.00000 1266.528 25.660004 44
## 59 59 2 60.30459 22.126696 142 128.3633 31.00000 1186.056 16.968849 44
## 60 60 2 34.73252 144.096893 141 121.0000 35.00000 1163.673 16.732523 44
## 61 61 2 35.36098 143.841354 140 120.3669 35.00587 1158.460 17.531569 44
## 62 62 2 12.00000 283.724396 140 112.0177 40.00000 1159.308 18.017662 44
## 63 63 2 28.63475 241.912872 140 115.0027 36.00000 1153.634 13.628904 44
## 64 64 2 38.00000 111.528107 141 123.0000 34.00000 1202.460 3.632200 44
## 65 65 2 19.92366 208.842209 141 117.0000 35.00000 1243.084 10.923658 44
## 66 66 2 16.00000 226.955856 140 115.9853 38.00000 1169.255 21.000000 44
## 67 67 2 38.26402 116.910995 141 122.9640 33.00000 1208.814 18.080193 44
## 68 68 2 53.00000 154.381088 139 120.0000 33.00000 1176.578 3.332412 44
## 69 69 2 47.64957 95.506989 141 124.0000 33.00000 1165.147 11.000000 44
## 70 70 2 34.96188 21.545456 142 128.0000 33.00000 1250.021 20.615139 44
## 71 71 2 17.35966 179.844604 140 118.0000 38.00000 1172.090 13.047443 44
## 72 72 2 30.00000 146.185760 140 120.0000 36.00000 1162.991 2.000000 44
## 73 73 2 34.34344 126.360092 141 122.0000 35.00000 1188.610 16.956320 44
## 74 74 2 38.02332 106.306831 141 123.0233 35.00000 1180.208 6.962867 44
## 75 75 2 54.00000 69.539719 142 126.0000 32.00000 1168.274 7.594716 44
## 76 76 2 26.98126 213.554977 141 117.0000 35.00000 1144.254 7.608913 44
## 77 77 2 82.23924 33.021381 140 127.9518 30.00000 1193.644 16.364378 44
## 78 78 2 37.62527 149.359497 141 120.0293 34.00000 1158.269 6.029291 44
## 79 79 2 53.95367 40.422386 141 127.9537 32.00000 1189.704 8.982679 44
## 80 80 2 66.00000 19.230110 142 129.0000 30.00000 1176.920 1.000000 44
## 81 81 2 17.00000 210.726776 140 117.0000 38.00000 1167.195 10.244825 44
## 82 82 2 10.00000 305.004395 140 112.0000 38.97747 1209.088 23.085899 44
## 83 83 2 71.25963 32.688908 142 128.0432 29.00000 1168.091 18.400412 44
## 84 84 2 16.59008 210.250671 140 117.0000 38.00000 1189.317 15.659211 44
## 85 85 2 18.44054 202.489563 140 117.0441 37.00000 1190.649 11.088174 44
## 86 86 2 48.94530 72.169975 142 126.0000 32.00000 1199.115 13.619845 44
## 87 87 2 36.69690 121.016777 141 121.3685 34.02003 1167.205 20.599169 44
## 88 88 2 63.77258 40.579601 142 128.0000 31.00000 1166.070 13.613713 44
## 89 89 2 35.73909 107.742851 141 123.0000 35.00000 1181.229 15.319310 44
## 90 90 2 29.37322 164.632294 141 120.0000 35.00000 1225.180 12.000000 44
## 91 91 2 44.42871 143.166260 141 121.0090 33.00000 1145.764 10.428706 44
## 92 92 2 66.23499 84.663139 142 125.0174 30.00000 1132.459 10.431952 44
## 93 93 2 20.62150 195.659027 141 118.0000 35.00000 1248.306 11.000000 44
## 94 94 2 16.00000 194.681168 140 117.9908 38.00000 1170.047 10.103127 44
## 95 95 2 42.03058 123.655266 141 122.0306 34.00000 1163.380 7.030577 44
## 96 96 2 40.01698 103.714996 141 124.0000 34.00000 1194.018 8.586416 44
## 97 97 2 65.00981 62.735703 142 126.0000 30.00000 1144.148 11.147452 44
## 98 98 2 49.03423 85.221207 141 125.0000 33.00000 1166.159 8.986865 44
## 99 99 2 79.17804 66.546135 142 126.0000 28.00000 1112.752 10.821962 45
## 100 100 2 27.96902 170.000351 141 119.0000 36.00000 1166.285 12.667546 44
## 101 101 2 21.32281 167.763794 140 119.0000 37.00000 1172.219 15.029562 44
## 102 102 2 54.71941 52.942146 142 126.3812 32.00000 1181.205 21.899353 44
## 103 103 2 33.38015 185.822006 141 118.0000 34.01198 1143.298 13.607867 44
## 104 104 2 28.37443 150.247620 141 120.0000 36.00000 1171.247 4.374433 44
## 105 105 2 30.33723 154.650818 141 120.0000 35.00000 1181.551 17.973600 44
## 106 106 2 50.09616 61.259510 141 126.0481 33.00000 1191.023 18.797804 44
## 107 107 2 39.42589 115.959381 141 123.0000 34.00000 1187.744 11.260800 44
## 108 108 2 21.90289 135.228958 141 121.9029 34.00000 1252.930 15.555657 44
## 109 109 2 14.01309 252.392471 140 115.0000 39.00000 1182.723 24.097075 44
## 110 110 2 44.95111 82.029579 141 125.0000 33.00000 1193.167 12.902219 44
## 111 111 2 14.00000 272.352417 139 114.0000 39.00000 1172.712 15.649142 44
## 112 112 2 86.42541 59.828053 142 127.0000 27.00000 1102.239 12.013014 45
## 113 113 2 37.62217 136.897842 141 121.0000 35.00000 1159.452 5.746593 44
## 114 114 2 54.73362 130.293320 141 122.0014 33.00000 1153.900 27.268824 44
## 115 115 2 41.23255 89.895660 142 124.9531 32.00000 1253.917 23.646996 44
## 116 116 2 79.00000 24.149178 142 129.0000 29.00000 1143.306 2.988217 44
## 117 117 2 38.36826 129.463318 141 122.0000 33.00000 1225.491 9.610627 44
## 118 118 2 34.37909 147.659988 141 121.0000 34.00000 1189.834 8.758190 44
## 119 119 2 78.59055 60.026829 142 126.0154 28.00000 1117.547 9.538243 45
## 120 120 2 47.77022 74.713654 141 125.6732 32.95148 1191.203 23.423868 44
## tseas vpd
## 1 55 35.00000
## 2 47 29.00000
## 3 57 31.00000
## 4 47 25.00000
## 5 45 28.00000
## 6 60 36.00000
## 7 48 31.00000
## 8 50 28.00000
## 9 55 31.00000
## 10 47 28.00000
## 11 56 31.00000
## 12 56 30.00000
## 13 51 30.00000
## 14 59 37.00000
## 15 58 34.00000
## 16 55 33.00000
## 17 46 26.00000
## 18 49 30.00000
## 19 51 32.04100
## 20 46 27.00000
## 21 53 30.00000
## 22 48 29.00000
## 23 58 32.00000
## 24 49 30.00000
## 25 52 30.00000
## 26 54 30.00000
## 27 52 30.00000
## 28 56 30.00000
## 29 60 34.00000
## 30 52 30.00000
## 31 52 29.00000
## 32 60 34.00000
## 33 58 32.00000
## 34 53 30.00000
## 35 56 34.00000
## 36 59 34.00000
## 37 51 29.00000
## 38 47 29.00000
## 39 NaN NaN
## 40 57 32.00000
## 41 58 32.00000
## 42 57 36.00000
## 43 54 33.00000
## 44 45 28.00000
## 45 46 26.00000
## 46 54 31.00000
## 47 57 32.00000
## 48 50 31.00000
## 49 59 29.00000
## 50 45 28.00000
## 51 46 29.00000
## 52 54 31.00000
## 53 48 29.00000
## 54 45 29.00000
## 55 55 31.00000
## 56 55 31.00000
## 57 45 29.00000
## 58 59 34.94745
## 59 60 33.00000
## 60 48 29.00000
## 61 46 28.00000
## 62 46 25.00000
## 63 46 26.00000
## 64 52 29.00000
## 65 46 29.00000
## 66 46 26.00000
## 67 53 31.00000
## 68 49 26.00000
## 69 54 31.00000
## 70 57 34.93906
## 71 45 27.00000
## 72 45 28.00000
## 73 51 29.00000
## 74 53 30.00000
## 75 56 32.00000
## 76 45 28.00000
## 77 59 29.00000
## 78 47 30.00000
## 79 59 32.00000
## 80 60 33.00000
## 81 46 26.00000
## 82 45 25.00000
## 83 59 34.00000
## 84 46 27.00000
## 85 46 27.00000
## 86 57 32.00000
## 87 52 30.00000
## 88 59 33.00000
## 89 52 29.00000
## 90 48 29.00000
## 91 48 30.00000
## 92 55 33.00000
## 93 46 29.00000
## 94 45 27.00000
## 95 52 30.00000
## 96 53 29.00000
## 97 56 32.00000
## 98 55 31.00000
## 99 56 34.00000
## 100 45 29.00000
## 101 45 28.00000
## 102 58 32.00000
## 103 45 29.00000
## 104 46 29.00000
## 105 48 29.00000
## 106 57 31.00000
## 107 53 30.00000
## 108 49 32.00000
## 109 47 26.00000
## 110 56 31.00000
## 111 47 25.00000
## 112 57 34.00000
## 113 47 28.00000
## 114 49 28.00000
## 115 55 33.00000
## 116 60 34.00000
## 117 51 30.00000
## 118 50 30.00000
## 119 56 34.00000
## 120 56 31.00000This is very useful in statistical model building (which we will get to in a week or two.)
Make a raster layer from vector data
To convert vector data to a raster we convert first to SpatVector then to raster with rasterize
plants %>% as("SpatVector") %>%
rasterize(dem) %>%
tm_shape() +
tm_raster(col = "dem", palette = "Greys")## Warning in showSRID(SRS_string, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum New Zealand Geodetic Datum 2000 in Proj4
## definition